Old Dominion University Desearch Foundatio: 


NASA-CR-190902 


DEPARTMENT OF MECHANICAL ENGINEERING & MECHANICS 
COLLEGE OF ENGINEERING & TECHNOLOGY 
OLD DOMINION UNIVERSITY 
NORFOLK, VIRGINIA 23529 


INVESTIGATION OF ADVANCING FRONT METHOD FOR 
GENERATING UNSTRUCTURED GRID 


By 

A. M. Thomas, Graduate Research Assistant 


//!/ — -O'/L 

/J/3V7 


and 

S. N. Tiwari, Principal Investigator 


f. 




Progress Report 

For the period ended May 31, 1992 


Prepared for 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, Virginia 23665 


Under 

Research Grant NCC1-68 

Dr. Robert E. Smith Jr., Technical Monitor 

ACD-Computer Applications Branch 


June 1992 


ro 



>* 



o 


r- 

4“ 


4* 

m 

fU 

n 

1 

f— 


f\l 

u 

fM 

O 

c 

rH 

Z 

D 

o 


<* 

n 

ro 

O 


u 

O r- 1 

ro 
z 

O C7^ 

t-H C ♦ 

K O — > 

< 

c cc cc c c 

^ a o v o 

IL 

uO Q13 C 

UJ C UJ o o 

> o o: ■* - 

Z X W C 

^ 1— h- <D •- 

UJ u Cc E 

r d o 
^ or * o 

nj h h 4J 

Q Z Irt U T 3 

o a z o 

O Of D O.C 

0 li ^ 

O a: 

1 o z 

z H cn M 

U ^ H ^ O 

i o < a a 

< Z OC W rH 

o < u. ct a 

< > Z O >* 
z a uj u rc 
v < o a z n 


DEPARTMENT OF MECHANICAL ENGINEERING & MECHANICS 
COLLEGE OF ENGINEERING & TECHNOLOGY 
OLD DOMINION UNIVERSITY 
NORFOLK, VIRGINIA 23529 


INVESTIGATION OF ADVANCING FRONT METHOD FOR 
GENERATING UNSTRUCTURED GRID 

By 

A. M. Thomas, Graduate Research Assistant 
and 

S. N. Tiwari, Principal Investigator 


Progress Report 

For the period ended May 31, 1992 


Prepared for 

National Aeronautics and Space Administration 
Langley Research Center 
Hampton, Virginia 23665 


Under 

Research Grant NCC1-68 

Dr. Robert E. Smith Jr., Technical Monitor 

ACD-Computer Applications Branch 



Submitted by the 

Old Dominion University Research Foundation 
P.O. Box 6369 

Norfolk, Virginia 23508-0369 


June 1992 


FOREWORD 


This is a progress report on the research project “ Numerical Solutions of Three-Dimensional 
Navier-Stokes Equations for Closed - Bluff Bodies” for the period ended May 31, 1992. Specific 
efforts during this period were directed in the area of “Investigation of Advancing Front Method 
for Generating Unstructured Grid.” 

This work was supported by the NASA Langley Research Center through Cooperative 
Agreement NCC1-68. The cooperative agreement was monitored by Dr. Robert E. Smith, Jr., 
of the Analysis and Computational Division (Computer Applications Branch), NASA Langley 
Research Center, Mail Stop 125. 
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SUMMARY 


Investigation of Advancing Front Method 
for Generating Unstructured Grid. 


A. M. Thomas* and S. N. Tiwari^ 

Old Dominion University, Norfolk, Virginia 23529 


In this report advancing front technique is used to generate an unstructured grid about simple 
aerodynamic geometries. Unstructured grids are generated using VGRID2D and VGRID3D 
software. Specific problem considered are NACA 0012 airfoil, bi-plane consisting of two NACA 
0012 airfoil, a four element airfoil in its landing configuration and an ONERA M6 wing. Inviscid 
time dependent solutions are computed on these geometries using USM3D and the results are 
compared with standard test results obtained by other investigators. A grid convergence study 
is conducted for NACA 0012 airfoil and compared with structured grid. Structured grid is 
generated using GRIDGEN software and inviscid solutions computed using CFL3D flow solver. 

The results obtained by unstructured grid for NACA 0012 airfoil showed an asymmetric 
distribution of flow quantities and a fine distribution of grid was required to remove this 
asymmetry. On the other hand structured grid predicted a very symmetric distribution but when 
the total number of points were compared to obtain the same results it was seen that structured 
grid required more number of grid points. 


Graduate Research Assistant. Department of Mechanical Engineering and Mechanics. 
Eminent Professor. Department of Mechanical Engineering and Mechanics. 
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1. INTRODUCTION 


A problem of great concern in computational fluid dynamics is the generation of suitable 
grids. Despite considerable effort devoted t oward the development for various grid generation 
methods, the process of developing a quality grid around a complex configuration and getting 
a flow solution still remains a challenging task, 

A grid is defined as a set of points with appropriate connections between the points. The 
points act as reference positions within the field at which the flow variables are to be computed 
and the connections between the points act as pathways for transferring information around the 
computational domain. In a structured grid the connectivity between the points is explicidy 
defined through a curvilinear coordinate system. In an unstructured grid the connectivity can 
be arbitrary and therefore must be specified. Solution methods that utilize a structured grid are 
generally more efficient than an unstructured grid solution method. However unstructured grids 
provide a much greater degree of flexibility than is available with a structured grid. In particular, 
unstructured grids can discretize a highly complex domain easily and are suitable for performing 
localized grid enrichment in solution adaptive methods. 

There are a variety of methods for generating unstructured grids. Among the different 
techniques are Delauny trianguiation methods [I], the modified Octree method (2) and the 
advancing front technique [3], In the present study, an advancing front technique is used to 
generate the unstructured grid for a variety of geometries. This method was selected because it 
does not require a separate library of modules to distribute grid points throughout the domain in 
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advance like the Delauny triangulation methods. The main disadvantage of the modified Octree 
method is that it has the difficulty of correctly defining the boundary segments. 

There are a number of software available to generate structured grid. Among the various 
codes available, GRIDGEN [4] code is used to generate structured grid around NACA 0012 
airfoil. This code provides an efficient means of developing multiple block grid for complex 
configurations. 

Steady state Euler solutions are computed using USM3D [5] on unstructured grids and 

CFL3D [6] on structured grids. The results obtained for NACA 0012 airfoil are plotted and 
compared. 

The objective of this report is to describe our experience with the advancing front algorithm 
for several test problems and to conduct a comparative study with structured grids. Steady state 
Euler solutions computed on these grids using the USM3D and CFL3D flow solver are presented. 
The grid generation technique is described in Sec. 2 and the flow solution methods are explained 
in Sec. 3. The numerical results are presented and discussed in Sec. 4. 
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2. GRID GENERATION 

2.1 Unstructured Grid Generation 


In this study unstructured grids are generated around simple 2-D and 3-D geometries. A 
2-D version of the unstructured grid generation package (VGRID2D) is used to generate the 
2-D grids. The input file contains the points defining the surface of the airfoil with connectivity 
table. A background grid is then setup which in 2-D is in the form of triangles. The background 
grid defines the local grid characteristic such as grid spacing. A typical background grid for a 
NACA 0012 airfoil is shown in Fig. 2.1. The 2-D grid is generated by first placing points on 
the boundary segments that define the solution domain, and then discretizing the interior region. 
A program has been developed (GEN3D) that uses the generated 2-D grid and duplicates it in 
the third direction, to form prism elements. Using the prism centroid, each prism is divided 
into eight tetrahedral cells. 

2.2 Structured Grid Generation. 

Structured grid is generated around NACA 0012 airfoil using GRIDGEN software. This soft- 
ware essentially consists of three main codes, GRIDBLOCK, GRIDGEN2D, and GRIDGEN3D. 
These are discussed below. 

GRIDBLOCK : This is used to develop blocking structures and set interblock connections. 
The surface definition for NACA 0012 airfoil is read in, which is the same as the one used for 
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the unstructured grid. A single block is used to define the whole domain. The 3-D lines are 
drawn which define the edges of the block. Computational directions and dimensions are added 

to the system. After the blocking strategy, flow boundary conditions are applied on all the faces 
of the block. 

GRIDGEN2D : This is used to generate surface grids on the six faces of the block. The grid 
is generated in components using the block edges. The block edges are then used as boundary 
condition for block faces and the resulting block faces are used as boundary condition for the 
block interior. Each of the four edges of a subface is generated before proceeding to the interior. 
The grid point spacing is defined when distributing points on these edges. The spacing around 
the airfoil is maintained approximately the same as that in the unstructured grid. Grid points 

interior to the subface are generated using the elliptic solver to achieve orthogonality near the 
surface of the airfoil. 

GRIDGEN3D : This program is used to distribute points within the interior of the block 
to generate volume grid. This is accomplished with the batch procedure of GRIDGEN3D and 
run on the Cray-super computer. 


3. GOVERNING EQUATIONS 
AND METHOD OF SOLUTION 


3.1 Structured and Unstructured Grids 

The governing equations are the 3-D unsteady Euler equations for inviscid compressible 
flow. For a bounded domain Q with a boundary dtl, the time dependent Euler equations in 
integral form can be written as 

y / j Q dv + j J F(Q) • hdS = 0 (3.1) 

n au 


r p 1 


' [pu, pv, pw) 

pu 


(pu 2 + P,puv,puw) 

< pv 

> and F(Q) = < 

(pvu,pv 2 -+ P,pvw ) 

pw 


(pwu, pwv, pw 2 + P) 

s p £ 0 > 


, (( pe 0 + P)u,(pe 0 + P)v,(pe 0 + P)w) , 


In the preceeding equations, p is the den sity, u, v and w are the x, y, z components of the 
velocity, e 0 is the total energy per unit volume and P is the pressure. The equations are non- 
dimensionalized by a reference density p^ and speed of sound a^. Assuming an ideal gas, 
the pressure is written as 


P = (7 - 1) ^e 0 - t;p(u 2 + v 2 + w 2 )^ 

where 7 represents the ratio of specific heats. 


(3.2) 
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The inviscid flow field is computed with a three-dimensional upwind flow solver (USM3D) 
developed at NASA/LaRC [4].The spatial discretization is accomplished with a cell centered 
finite volume formulation using the flux difference splitting procedure. The solution is advanced 
in time using a 3-stage Runge-Kutta time stepping scheme. Local time stepping and implicit 
residual smoothing are used to accelerate the convergence of the solution to a steady state. 

For the structured grid the time dependent Euler equations are similar to Eq. (3.1). The 
unsteady three dimensional equations are solved by a 3-factor implicit time advancement 
algorithm (CFL3D) which uses a block tridiagonal inversion. Here also local time stepping 
and implicit residual smoothing are used to accelerate the convergence of the solution to a 
steady state. 

Boundary Conditions 

For the solid boundaries the flow tangency condition is imposed. A condition of zero mass 
and energy flux through the surface is ensured. Characteristic boundary conditions are applied to 
the far field subsonic boundary. At an outflow boundary, the two tangential velocity components 

and the entropy are extrapolated from the interior while at an inflow they are specified as having 
far- field values. 
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4. RESULTS AND DISCUSSION 


Results obtained for different cases and conditions are presented in this section. 

Case 1 

A NACA 0012 airfoil has been used to validate the algorithm developed. The program 
considers the 2-D grid points and extends them to a 3-D configuration compatible to the flow 
solver. The computational domain is bounded by a rectangular box with boundaries at -10 £ 
X <; 10, 0 <; y <; 2.0, -10 ^ z ^ 10, where the airfoil chord length is one. Initially a coarse 
grid was generated from VGRID2D consisting of 112 cells, 585 points with 53 boundary points 
(Fig. 4.0). After going through the 3-D algorithm in which five planes were specified, the 3-D 
configuration had 35584 cells, 7373 points with 1344 boundary points. 

The computations were performed at a CFL number of 4.0 and Mach number of 0.5 with 
zero angle of attack. For this case, the flow field should be symmetrical about the airfoil and 
no shock should occur. The solutions were started from free stream initial conditions with the 
first order scheme. When the RMS average of all residuals dropped by one order of magnitude 
the solver was switched on to the higher order scheme. 

Figure 4. 1 shows the convergence history with a decrease of approximately 2.5 orders of 

magnitude. This figure also shows the total lift and drag which are seen to be constant throughout 
the computational period. 

The percentage of total pressure loss over the surface is shown in Fig. 4.2. It can be 
seen that the pressure losses are symmetrical about the airfoil. The maximum magnitude of the 
pressure loss was observed to be 0.04 %. Theoretically, there should be no loss of total pressure 
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Fig. 4.1 Convergence history for the coarse grid. 








for this test case because the flow is an isentropic process. The observed pressure losses are 
due to numerics alone. 

Figures 4.3 show the Mach number contours. It can be seen from the figure that the 
computed flow field is symmetric and free of shocks. 

Figure 4.4 shows the coefficient of pressure on the upper and lower surfaces. Since there was 
no shock, the two curves are approximately symmetrical and no abrupt changes are observed. 
Theoretically the two curves should lie on top of each other because of the symmetric flow 
conditions. The deviation in the curve may be due to the coarseness of the grid. 

To investigate the asymmetric Cp distribution, a parametric study is conducted and compared 
with the results from the structured grid. A medium fine grid is selected which has 76 points on 
the airfoil and a total of 1583 points with 3054 cells. Figure 4.5 shows the plot of the grid in 
2-D. The initial condition for this test case is obtained by linearly interpolating the converged 
coarse grid solution onto the new grid. Figure 4.6 shows the convergence history for this case. 
The solution was started with a second order differencing scheme but after 150 iterations the 
solution appeared to be diverging. At this point the solver was switched back to the first order 
scheme and the solution converged with approximately three orders of magnitude reduction in the 
residuals. Figure 4.7 shows the Mach number contour plot. Figure 4.8 shows the Cp variation 
over the upper and lower surfaces. In the Cp plot there is now much better agreement between 
the upper and lower surface results. 

Next a very fine grid is selected in which 163 points are located on the airfoil and the 2-D 
grid contains a total of 2582 points and 4961 cells. Figure 4.9 shows the grid distribution. The 
converged solution from the medium fine grid is interpolated on to the fine grid to provide a 
better initial condition for this test case. In accordance with the approach used in the previous 
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Fig. 4J5 Medium fine mesh for NACA 0012 airfoil. 
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test case, the solution was obtained with a first order scheme. Figure 4.10 shows the convergence 
history for the entire parametric study. Figure 4. 1 1 shows the Cp distribution. For this test, the 
curves in the Cp plot have very good agreem ent. The Mach number contour plot is similar to 
the medium fine grid test case and hence not shown here. 

To obtain a good resolution of the contour plots, the grids were concentrated near the nose 
and tail of the airfoil. The 2-D grid contained 120 points on the airfoil and a total of 2193 points 
and 4229 cells. Figure 4.12 shows the grid distribution. The convergence history is illustrated in 
Fig. 4.13. The solution was started from a free stream initial condition with a first order scheme, 
and when the residuals dropped by one order of magnitude the solver automatically switched over 
to a second order scheme. A steady state solution was obtained with approximately two orders 
of magnitude reduction in the residuals. Figure 4.14 shows the Mach number plot where a clear 
resolution is seen near the nose and tail of the airfoil. Figure 4.15 shows the Cp distribution. 
Comparison of pressure distribution for the four different cases is shown in Figs. 4.16 and 4.17. 
The results show that the solution improves with the increase in the number of grid points on 
the airfoil. Hence it is concluded that a very close resolution of the grid point is required near 
the airfoil to attain a satisfactory result. 

CASE 2 

To confirm the results obtained from the unstructured grid, and to have a comparative study, 
structured grid was generated around the same airfoil with GRIDGEN code. 

The physical dimensions were maintained the same as that in the unstructured grid. A C-H 
type grid was generated and computations were performed with a Mach number of 0.5 at zero 
angle of attack. 
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Fig. 4.11 Variation of coefficient of pressure over the upper and lower surface. 
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Fig. 4.12 Adapted grid for NACA 0012 airfoil. 
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contours for the adapted grid. 
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Fig. 4.15 Variation of coefficient of pressure over the upper and lower surface. 
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Fig. 4.16 Comparison of Cp on the upper surface for four different grid concentration. 
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Fig. 4.17 Comparison of Cp on the lower surface for four different grid concentration 
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Fig. 4.18 Initial coarse structured grid for NACA 0012 airfoil. 



Figure 4.18 shows the coarse grid with 38 points around the airfoil and 49 points in the 
normal direction. Figure 4.19 shows the Mach number contour plot and Fig. 4.20 the Cp 
distribution. Both figures show the symmetric distribution of the flow field quantities. The 
Cp values for the top and bottom surfaces exactly match each other. Figure 4.21 shows the 
comparison of Cp between structured and unstructured grids. Near the stagnation point, the 
unstructured grid predicts the flow accurately but the structured grid gives a poor result. This 
might be due to the coarse distribution of the grid points near the leading edge. 

To investigate the leading edge inaccuracy, a medium fine grid is selected with the same 
number of grid points as in the unstructured grid i.e., 76 points. Figure 4.22 shows the grid 
distribution. Computations were performed with free stream initial conditions, and Mach number 
contours for a steady state solution are plotted in Fig. 4.23. Again, a symmetric distribution 
of flow field quantities is seen. The Cp plot is shown in Fig. 4.24 and Fig. 4.25 shows the 
comparison with the unstructured grid. A better prediction of the flow field quantities is seen 
near the leading edge. 

Next, a very fine grid with 163 points around the surface of the airfoil is generated. In this 
case also the solution was started with the free stream initial condition. Figure 4.26 shows the 
grid distribution and Fig. 4.27 the Mach number contour plot. Coefficient of pressure for the 
top and bottom surface is shown in Fig. 4.28. The flow field is symmetrical and a very good 
estimation of the flow field quantities is seen near the leading edge. Figure 4.29 shows the Cp 
comparison with the unstructured grid and Fig. 4.30 shows the comparison of Cp for different 
grid concentration. It is seen that as the grid is concentrated around the airfoil better results 
are obtained. Coefficient of pressure distribution along the surface for both structured and 
unstructured grid is compared with the experimental values obtained by Harris [9] in Fig. 4.31. 
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Fig. 4.19 Mach number contours for the coarse grid. 
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431 Comparison of Cp between lower and upper surface. 
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Fig. 4.22 Medium fine grid for NACA 0012 airfoil. 
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Fig. 4.23 Mach number contours for medium grid. 
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Fig. 4.25 Comparison of Cp between structured and unstructured grid. 
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Fig. 4.26 Fine grid for NACA 0012 airfoil. 




Fig. 4.27 Mach number contours for fine grid. 











Its seen that the structured grid has a much better agreement with the experimental values than 
unstructured grid. This might be due to the efficient flow solver used for structured grid than 
the unstructured grid. 

Comparison between the structured and unstructured grid for the various grid concentration 
is shown in Table 1. The table also shows the percentage of error in pressure calculated at the 
stagnation point. It is seen from the table that structured grid requires more number of points to 
accurately predict the flow field when compared to unstructured grid. 

Case 3 

In this case the same airfoil was considered as in case 1 but for a Mach number of 0.84 and 
1.25° angle of attack. Figures 4.32 and 4.33 show the Mach number and pressure contours. Here, 
a shock in the form of an expansion wave is observed over the top surface and a compression 
wave at the bottom surface. Figure 4.34 shows the coefficient of pressure plotted over the upper 
and lower surface. For this case the lift and drag coefficient were computed to be Q = 0.07, 
C d = 0.02 
Case 4 

In this test case a biplane is considered made up of two NACA 0012 airfoils. The outer 
boundaries are located at a distance of ten chord lengths. The distance between the two airfoils 
is one chord length with the top airfoil shifted from the center by 0.5 chord lengths. The grid 
for this case is shown in fig. 4.35. The above geometry was tested for a Mach Number of 0.8 
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Fig. 4.29 Comparison of Cp between structured and unstructured grid for a fine mesh. 
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Fig. 431 Comparison of Cp with experimental values. 
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Fig. 4.35 Unstructured grid for a bi-plane. 













and zero angle of attack. Overall, the Mach number contour plot in Fig. 4.36 is in good 
agreement with the results obtained by Hwang [7]. 

Case 5 

In this case a four-element airfoil system is considered. The geometry represents a section 
of a wing in its landing configuration. Figure 4.37 shows the coarse mesh generated in 2-D with 
5354 points, 10469 cells and 245 boundary points. Figure 4.38 shows the fine mesh consisting 
of 13898 points, 454 boundary points and 27348 cells. 

The coarse mesh was taken for parametric study and run under the same conditions as that 
in Ref. [1] ( ie. a Mach Number of 0.3 and 5 degrees angle of attack ). Figure 4.39 shows the 
convergence history with first order differencing. The residuals dropped by about one and half 
orders of magnitude and then remained constant. Figure 4.40 shows the convergence history 
when the program was set for automatic selection of differencing. The program started with first 
order differencing and when the residuals dropped by one order of magnitude it switched over to 
second order. The jump in the curve shows the change in the order of differencing. The second 
order differencing however did not work as expected since a very coarse mesh was selected. 
The plots shown for the variation of different parameters are plotted for first order differencing. 

Figure 4.41 shows the coefficient of pressure plotted over the top and bottom surfaces. The 
figure has the same trend as that shown in Ref. [1], The variation in values are due to the 
fact that in Ref. [1] solution adaptive techniques are employed and hence a better grid point 
distribution is obtained. 

Figures 4.42 and 4.43 show the Mach number and pressure contours. The plots are in good 
agreement with that shown in Ref. [8], 
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Fig. 437 Initial coarse mesh for a 4-element airfoil in its landing configuration. 
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Case 6 


An ONERA M6 wing is used as a 3-D test case to generate unstructured grid. The wing has 
a leading edge sweep of 30 degrees, an aspect ratio of 3.8, taper ratio of 0.56, and symmetrical 
airfoil sections. The wing has a root chord of 0.67 and a semispan b of 1.0 with a rounded 
tip. The computational domain is bounded by a rectangular box with boundaries at - 6.5 <; 
x £ 11.0, 0.0 £ y <; 2.5 and - 6.5 <, z <, 6.5. Transonic solutions are computed on the two 
grids shown in Figs. 4.44 and 4.45 at Moo = 0.84 and a = 3.06°. The solutions were started 
with a free stream initial conditions with first order scheme and when the residuals dropped 
by one order of magnitude, at which the solver automatically switched over to higher order 
scheme. A comparison of wing surface pressure contours for the two mesh is shown in Figs. 
4.46 and 4.47. The contours on both meshes show double shock wave on the upper surface. 

The contours on mesh2 show much sharper shock waves than the first one, demonstrating an 
effect of mesh density. 

Coefficients of lift drag and pitching moment about the wing apex are given in Table 2. The 

table also shows the number of grid points and cells in each of the mesh used. The values are 
compared with that in Ref. [10], 
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Fig. 4.44 Initial coarse mesh for ONER A M6 wing. 
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Fig. 4,45 Fine mesh for ONERA M6 wing- 
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Mesh 1 

Mesh 2 

Number of Cells 

45,335 / 47,344* 

94,257 / 98,317* 

Number of Nodes 

9,062 / 9,401* 

18,357 / 19,048* 

Number of Boundary Faces 

2,829 / 2,932* 

5,146 / 5,196* 

Lift Coefficient 

0.2879 / 0.2892* 

0.2912 / 0.2893* 

Drag Coefficient 

0.0184 / 0.0195* 

0.0154 / 0.0167* 

Moment Coefficient 

-0.1726 / -0.1715* 

-0.1733 / -0.1702* 


* Values taken from Ref. 

Table 2. Comparison of ONERA M6 wing values with Ref. [10] 
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5. CONCLUDING REMARKS 


In this study the superiority of unstructured grid over structured grid is demonstrated. 
Unstructured grids generated by VGRID3D software is compared with structured grids generated 
by GRIDGEN code. Specific problems considered are NACA 0012 airfoil, bi-plane consisting 
of two NACA 0012 airfoils, a four element airfoil in its landing configuration and an ONER A 
M6 wing. Inviscid time dependent solutions are computed for these geometries and compared 
with standard test results. A grid convergence study is conducted for NACA 0012 airfoil and 
compared with structured grid. The results obtained by NACA 0012 airfoil using unstructured 
grid showed an asymmetric distribution of flow field quantities and a fine distribution of grid 
was required to remove this asymmetry. On the other hand, the structured grid predicted a very 
symmetric distribution, but when the total number of points were compared the structured grid 
required more number of points. 

For future studies its is recommended that unstructured grids be used with improved surface 
definition like (NURBS) around complex configurations. It is also advocated that refined model 
geometry be created before analyses are performed. Three-dimensional prismatic elements which 
are structured in the normal direction and unstructured in the axial direction will be a good test 
case for incorporating viscous terms. 
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